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Abstract.  We  consider  a  discrete  time-and-space  route-optimization  problem,  across  a  fi¬ 
nite  time  horizon,  in  which  multiple  searchers  seek  to  detect  one  or  more  probabilistically 
moving  targets.  The  paper  formulates  a  novel  convex  mixed-integer  nonlinear  program  for 
this  problem  that  generalizes  earlier  models  to  situations  with  multiple  targets,  searcher 
deconfliction,  and  target-  and  location-dependent  search  effectiveness.  We  present  two  so¬ 
lution  approaches,  one  based  on  the  cutting-plane  method  and  the  other  on  linearization. 
These  approaches  result  in  the  hrst  practical,  exact  algorithms  for  solving  this  important 
problem,  which  arises  broadly  in  military,  rescue,  law  enforcement,  and  border  patrol  op¬ 
erations.  The  cutting-plane  approach  solves  many  realistically  sized  problem  instances  in 
few  minutes,  while  existing  branch-and-bound  algorithms  fail.  A  specialized  cut  improves 
solution  times  by  50%  in  difficult  problem  instances.  The  approach  based  on  linearization, 
which  is  applicable  in  important  special  cases,  may  further  reduce  solution  times  with  one  or 
two  orders  of  magnitude.  The  solution  times  for  the  cutting-plane  approach  tend  to  remain 
constant  as  the  number  of  searchers  grows.  In  part,  then,  we  overcome  the  difficulty  that 
earlier  solution  methods  have  with  many  searchers. 
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1  Introduction 


We  consider  a  discrete-time  route-optimization  problem,  denoted  SP,  in  which  multiple, 
fully  cooperating  searchers  and  one  or  more  probabilistically  moving  targets  operate  in  an 
environment  consisting  of  a  hnite  set  of  cells.  In  the  case  of  a  single  target,  the  searchers 
seek  to  minimize  the  probability  of  not  detecting  the  target  across  a  time  horizon.  In  the 
presence  of  multiple  targets,  the  searchers  adopt  a  worst-case  approach  and  aim  to  minimize 
the  probability  of  not  detecting  the  target  with  the  largest  nondetection  probability.  The 
searchers  are  subject  to  constraints  on  route  continuity,  may  require  some  travel  time  to 
move  from  one  cell  to  another,  and  may  need  to  maintain  a  minimum  distance  between  each 
other  to  avoid  interference.  The  searchers  are  imperfect  and  may  not  detect  a  target  present 
in  a  cell  subject  to  search.  However,  we  ignore  the  possibility  of  false  positive  reports  of 
targets. 

SP  arises  broadly  in  military,  rescue,  law  enforcement,  and  border  patrol  operations.  In 
military  operations,  the  searchers  may  be  manned  or  unmanned  aircraft  looking  for  suspected 
individuals  or  downed  pilots  in  an  area  of  interest.  Park  rangers  may  search  for  lost  hikers. 
In  a  damaged  or  burning  building,  hre  hghters  and  ground  robots  may  search  for  trapped 
individuals.  Law  enforcement  officers  may  act  as  searchers  when  looking  for  criminals.  Near 
national  borders,  the  searchers  are  border  patrols  seeking  illegal  immigrants.  The  searchers 
may  also  be  Coast  Guard  cutters  and  helicopters  scanning  the  ocean  for  ship  wrecks  and 
smugglers.  We  adopt  the  traditional  terminology  (see,  e.g.,  [33,  31])  and  refer  to  the  objects 
being  looked  for  as  “targets,”  even  though  the  searchers  may  not  wish  those  objects  any 
harm.  Often  in  these  applications,  the  targets  are  unaware  of  the  searchers’  routes  or  are 
unable  to  take  advantage  of  such  information.  Hence,  it  is  realistic  to  assume  that  the  targets 
move  between  cells  according  to  some  known  probability  distribution  that  is  independent 
of  the  searchers’  routes;  see  Section  2  for  details.  This  assumption  is  also  found  in,  e.g., 
[29,  9,  33,  17]  and  we  adopt  it  throughout  this  paper.  The  possibility  of  a  target  intelligently 
adapting  to  the  searchers’  routes  is  beyond  the  scope  of  the  paper. 
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The  situations  described  above  may  involve  just  a  few  searchers.  However,  with  the  tech¬ 
nological  improvement  of  autonomous  systems  and  sensors,  planners  now  frequently  need  to 
optimally  route  many  searchers.  (We  refer  to  [25]  for  an  overview  of  recent  advances  in  the 
area  of  autonomous  aerial  systems  and  [20]  for  an  introduction  to  the  U.S.  Department  of 
Defense’s  effort  to  utilize  a  large  number  of  autonomous  systems  for  reconnaissance,  infor¬ 
mation  gathering,  and  search  and  rescue  operations.)  Effective  manual  planning  for  a  few 
hnman-controlled  searchers  is  difficulty  (see  for  example  p.7-7  in  [33])  and,  of  conrse,  even 
more  difhcnlt  for  a  large  nnmber  of  autonomous  systems.  In  fact,  [20],  p.  164,  states  that 

“Many  recognize  the  potential  of  robotics/nnmanned  systems  in  military  opera¬ 
tions;  yet  few  nnderstand  how  to  employ  them  effectively.  Computer  generated 
modeling  and  simnlation  tools  are  extremely  valuable  in  an  attempt  to  generate 
both  operational  concepts  for  tactics,  techniqnes  and  procedures  (TTP)  as  well 
as  technical  requirements  to  enable  those  procedures.” 

Solutions  of  SP  recommend  effective  rontes  for  the  searchers  that  help  planners  to  better 
utilize  their  assets,  provide  input  to  control  units  of  antonomous  systems,  and  possibly  save 
lives,  for  example,  in  search  and  rescue  operations.  The  U.S  Coast  Guard’s  tactical  decision 
aid  GASP  and  the  U.S.  Navy’s  NODESTAR  both  ntilize  solutions  of  special  instances  of  SP 
obtained  by  henristics  (see  Chapter  7  in  [33]). 

Previous  studies  have  focused  on  single-searcher  and  single-target  instances  of  SP;  see 
[26]  for  a  comprehensive  review  of  papers  np  to  1991;  [17,  28],  and  references  therein,  cover 
more  recent  work.  The  most  efficient  exact  algorithms  for  such  instances  of  SP  appear  to  be 
specialized  branch-and-bonnd  algorithms.  These  algorithms  obtain  bounds  by  replacing,  in 
effect,  the  probability  of  detection  with  the  expected  nnmber  of  detections  [7,  32,  17,  28]. 

In  principle,  such  specialized  branch-and-bound  algorithms  are  also  applicable  to  in¬ 
stances  of  SP  with  mnltiple  searchers  [7,  27].  However,  when  applied  to  realistically  sized 
instances  with  more  than  two  searchers,  they  tend  to  reqnire  excess  compnter  memory  and 
prohibitive  solution  times  [27].  This  has  motivated  the  development  of  heuristics  for  SP, 
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and  related  problems,  such  as  local  search  and  genetic  algorithms  [7]  as  well  as  the  cross¬ 
entropy  method  [27].  Other  heuristic  algorithms  rely  on  approximations  of  probabilities 
by  expectations  [7],  myopic  optimization  with  a  receding-horizon  [7,  12,  35,  23],  sequential 
optimization  of  each  searcher  [29,  35,  13],  and  decentralized  optimization  by  each  searcher 
[3,  36,  35].  Of  these  studies,  only  [36,  35]  deal  with  multiple  targets,  and  then,  only  in  the 
context  of  heuristics.  Hence,  there  is  a  need  to  model  multi-searcher  instances  of  SP  and  to 
provide  exact  algorithms  for  their  solution. 

This  paper  formulates  a  novel  convex  mixed-integer  nonlinear  program  (MINLP)  for 
SP.  The  program  generalizes  the  formulation  in  [29]  to  multiple  targets.  It  also  appears  to 
be  the  hrst  to  prevent  interference  between  searchers  (i.e.,  to  “deconflict”  searchers)  in  this 
context  and  to  account  for  the  fact  that  a  searcher’s  effectiveness  may  depend  on  the  target, 
the  time  of  search,  and  the  searcher’s  previous  location.  The  latter  generalization  provides 
an  important  practical  advantage  as  it  reduces  the  need  for  hne  time  discretizations  when 
modeling  real-world  search  missions. 

We  present  two  solution  approaches  to  MINLP  resulting  in  the  hrst  practical,  exact  al¬ 
gorithms  for  SP  with  multiple  searchers.  One  approach  based  on  the  cutting-plane  method 
(see,  e.g.,  [15,  8,  34]),  leads  to  fast  solutions  for  many  realistically  sized  problem  instances 
for  which  existing  branch-and-bound  algorithms  fail.  For  difficult  instances  of  SP,  we  im¬ 
prove  solution  times  further  with  a  specialized  cut.  The  other  approach  is  based  on  novel, 
equivalent  linearizations  of  MINLP  available  in  important  special  cases.  Optimal  solutions 
of  the  linearizations  are  often  easily  obtained  by  standard  mixed-integer  linear  programming 
solvers  leading  to  further  reduction  in  solution  times  by  one  to  two  orders  of  magnitude. 

Poor  scalability  of  algorithms  as  the  number  of  searchers  grows  is  of  major  concern  when 
solving  SP  and  related  problems  and  has  led  to  the  development  of  heuristic  algorithms  (see, 
e.g.,  [35,  13,  27]).  The  cutting-plane  approach  scales  well  as  the  solution  times  tend  to  remain 
constant  as  the  number  of  searchers  grows.  We  therefore  overcome,  in  part,  the  difficulty  of 
solving  problem  instances  with  many  searchers. 

The  remainder  of  the  paper  is  outlined  as  follows.  The  next  section  dehnes  SP  precisely. 
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Section  3  derives  mathematical  programming  formulations  of  SP  and  its  linearizations.  In 
Section  4,  we  present  algorithms  and  numerical  results,  with  particular  focus  on  an  important 
special  case  involving  a  single  target. 

2  Problem  Statement 

2.1  Searchers,  Targets,  and  their  Environment 

We  let  the  area  of  interest  (AOl)  be  discretized  into  a  hnite  set  of  cells  C  =  {!,... ,C},  and 
let  the  time  horizon  be  discretized  into  a  hnite  set  of  time  periods  Tq  =  {0}  U  T,  where 
T  =  {1,2,  ...,T}.  Search  for  targets  takes  place  during  time  periods  t  E  T,  with  t  =  0 
representing  the  time  period  prior  to  start  of  the  search.  There  are  K  independent  targets 
present  in  the  AOl  with  each  target  /c  G  /C  =  (1,2, ...,  iP}  occupying  one  cell  in  each  time 
period.  The  quantity  Uk^t  G  C  denotes  the  (random)  cell  that  target  k  occupies  during  time 
period  t  E  T.  The  vector  of  cells  Uk  =  {ujk,i,ujk,2,  ■  ■  ■  denotes  a  possible  path  for 

target  k,  and  qk{oJk)  denotes  the  given  probability  that  target  k  takes  that  path.  The  set  Vtk 
denotes  the  collection  of  all  possible  paths  for  target  k  with  positive  probability  qk{oJk)-  In 
practice,  the  data  klk  and  quiojk)  are  generated  using  Monte  Carlo  sampling  from  (complex) 
target  motion  models  (as  in  U.S.  Coast  Guard’s  decision  aid  GASP)  or  dehned  implicitly 
by  Markov  transition  matrices  (as  in  the  case  of  U.S.  Navy’s  decision  aid  NODESTAR);  see 
Chapter  7  in  [33] .  This  paper  considers  both  ways  of  specifying  target  movement  and  refers 
to  the  former  way  as  a  conditional  target  model  and  to  the  latter  way  as  a  Markovian  target 
model.  In  Sections  3  and  4,  we  see  that  it  is  not  necessary  to  explicitly  enumerate  all  possible 
target  paths  in  the  case  of  a  Markovian  target  model. 

It  is  trivial  to  extend  the  current  framework  to  situations  with  targets  that  may  not 
be  present  in  the  AOl  and  targets  that  enter  and  leave  the  AOl  during  the  time  horizon. 
However,  we  do  not  examine  that  situation  further. 

There  are  L  classes  of  searchers,  with  each  class  I  E  C  =  (1,2,  ...,L}  containing  J; 
identical  searchers.  During  each  time  period  t  E  Tq,  each  searcher  occupies  a  cell  or  is  in 
transit  between  cells.  When  occupying  a  cell  c,  a  searcher  of  class  /  may  select  to  move  to 
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any  cell  “adjacent”  to  c  as  defined  by  the  forward  star  J^i{c)  C  C.  We  also  let  7li{c)  C  C 
denote  the  reverse  star  of  cell  c,  which  represents  the  set  of  cells  from  which  a  searcher  of 
class  I  can  reach  cell  c  in  one  move.  By  convention,  c  G  J^i{c)  and  c  G  7^/(c).  A  searcher  of 
class  /  requires  di^c,c'  time  periods  to  move  from  cell  c  to  cell  c'  G  J^i{c)  and  to  search  cell  c' 
for  one  time  period.  Since  the  time  to  search  the  “destination”  cell  c'  is  included  in  di^c,c', 
di,c,c'  >  1  for  all  I,  c,  d  and  di^c,c'  =  1  only  if  the  time  to  move  from  c  to  d  is  zero. 

Searchers  may  interfere  with  each  other  and  could  be  required  to  maintain  a  minimum 
internal  distance.  We  let  Uc  be  the  maximum  number  of  searchers  allowed  to  occupy  cell  c 
during  any  one  time  period  t  G  T.  Moreover,  for  each  possible  move  between  to  cells  for 
a  searcher,  we  define  a  corresponding  set  of  incompatible  moves  between  cells  that  would 
cause  interference  if  carried  out  by  another  searcher.  Specihcally,  if  a  searcher  of  class  /  moves 
from  cell  c  to  cell  d  starting  in  time  period  t,  then  the  set  V{l,c,d,t)  gives  all  quadruples 
of  searcher  classes,  cell  pairs,  and  time  periods  that  are  incompatible  with  that  searcher’s 
move.  We  refer  to  these  restrictions  as  deconfliction  constraints. 

We  let  Xi^c,c',t  denote  the  number  of  searchers  of  class  I  that  occupy  cell  c  in  time  period 
t  eTq  and  that  move  to  cell  d  next,  and  let  X  denote  the  vector  with  components  Xi^c,c',t, 
I  E  C,  c,  d  E  C ,  and  t  &  Tq.  We  refer  to  X  as  a  search  plan. 

2.2  Sensor  Model 

We  assume  that  each  searcher  is  equipped  with  one  imperfect  sensor.  Each  time  period  t  G  T 
in  which  a  searcher  occupies  a  cell,  the  searcher’s  sensor  takes  one  “look”  in  the  cell  for  each 
target.  When  a  searcher  is  in  transit  between  cells,  the  sensor  is  inactive.  The  probability 
that  one  look  for  a  target  in  a  cell  detects  the  target,  given  that  the  target  currently  occupies 
the  cell,  may  depend  on  the  searcher  class  (is  it  a  high-  or  low-quality  searcher?),  the  target’s 
characteristic  (is  it  shiny  or  camouflaged?),  the  cell  (is  it  forested  or  open?)  and  time  of  day 
(is  it  bright  mid-day  or  dark  midnight?).  Specifically,  if  target  k  and  a  searcher  of  class  I 
occupy  cell  c  in  time  period  t  and  d  is  the  searcher’s  previous  cell,  then  the  probability  that 
the  searcher’s  look  during  time  period  t  detects  the  target  is  gi^c',c,t,k  ^  (0, 1).  We  refer  to 
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this  probability  as  the  glimpse-detection  probability. 

We  note  that  the  glimpse-detection  probability  depends  on  the  searcher’s  previous  lo¬ 
cation.  This  dependence  may  arise  if  adjusting  search  pattern  and/or  altitude,  refocusing 
a  sensor,  and  becoming  familiar  with  a  new  cell  have  a  signihcant  detrimental  effect  on  the 
searcher’s  capability  to  detect  a  target.  In  addition,  this  dependance  allows  us  to  account 
indirectly  for  small  transit  times  (much  less  than  the  length  of  a  time  period)  between  cells 
by  reducing  the  glimpse-detection  probability  from  its  nominal  value  if  the  searcher  just 
moved  into  a  cell.  For  example,  suppose  that  the  real-world  travel  time  from  cell  d  to  c 
is  one  minute.  To  model  this  situation  (approximately),  we  would  normally  require  a  time 
period  of  (approximately)  one-minute  duration.  However,  this  may  result  in  a  large  number 
of  time  periods  and  long  computing  times.  Alternatively,  we  can  dehne  a  longer  time  period, 
say  10  minutes,  and  let  the  glimpse-detection  probability  in  cell  c  be  somewhat  reduced  if  a 
searcher’s  previous  cell  were  d  as  compared  to  if  it  were  c.  This  will  reflect  the  fact  that  a 
searcher  coming  from  cell  d  has  only  nine  minutes  to  search  c  compared  to  10  minutes  if  the 
searcher  had  already  been  present  in  c.  Hence,  we  avoid  adopting  a  hue  time  discretization 
with  resulting  high  computational  cost. 

We  assume  that  all  the  searchers’  looks  are  independent.  Hence,  given  search  plan  X, 
the  probability  that  no  searcher  detects  target  k  in  cell  c  in  time  period  f,  given  that  target 
k  occupies  cell  c  at  that  time,  equals 


l^C  c'g7?.;(c) 


I  ^  ^  '  ^l,c' I  ) 

V  l£Cc'&Tli{c) 

where  for  all  /  G  £,  c  G  C,  d  G  7^«(c),  t  G  T,  and  A;  G  /C, 


(1) 

(2) 


,c,t,k  ln(l  Ql^d 


(3) 


is  the  detection  rate.  While  the  glimpse-detection  probability  may  be  difficult  to  estimate 
directly  without  extensive  held  testing,  the  detection  rate  can  often  be  related  to  a  searcher’s 
speed  and  sensor  range,  the  size  of  the  cell,  and  the  length  of  the  time  period;  for  example 
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see  p.  2-1  in  [33]. 


SP  seeks  to  minimize,  by  choice  of  a  search  plan  X,  the  probability  of  not  detecting 
the  target  with  the  largest  nondetection  probability  during  the  time  horizon.  The  choice 
of  search  plan  is  subject  to  the  route  constraints  induced  by  the  forward  and  reverse  stars 
jF/(c)  and  Tliic),  deconfliction  constraints  given  by  Uc  and  V{l,c,c' ,t),  and  the  given  initial 
condition  that  xi^c,o  searchers  of  class  I  occupy  cell  c  in  time  period  0. 

3  Models  of  Search  Problem 

In  this  section,  we  formulate  SP,  in  complete  generality,  as  a  convex  MINLP.  Since  the 
model  is  nonlinear,  we  anticipate  relatively  long  solution  times  for  standard  MINLP  solvers 
(e.g.,  Bonmin  [6],  DICOPT  [11]).  Hence,  a  main  focus  of  this  paper  is  to  develop  solution 
approaches  that  utilize  structure  present  in  important  classes  of  problem  instances.  Conse¬ 
quently,  we  go  on  to  construct  two  linear  models  for  classes  of  problem  instances  involving 
homogenous  searchers  and  a  single  target. 

This  section  hrst  states  the  nonlinear  model  of  SP  and  second  deals  with  the  lineariza¬ 
tions  for  special  classes.  The  section  ends  with  a  discussion  of  linearization  of  the  full 
nonlinear  model. 

3.1  Nonlinear  Model  of  SP 

We  state  SP  as  a  convex  MINLP,  generalizing  the  formulation  in  [29]  to  account  for  multiple 
targets  and  deconfliction  constraints  as  well  as  glimpse-detection  probabilities  that  may 
depend  on  the  target,  the  time  of  search,  and  a  searcher’s  previous  location.  The  resulting 
program  takes  the  following  form. 

Model  SPX; 


Indices 

c,  c',  c",  d" 

cells  (c,  c',  c",  d"  G  C  =  {1, . . . ,  C}). 

t,  t' 

time  periods  {t,  t'  eTq  =  {0}  \JT 

Id' 

searcher  class  (/,  /'  G  £  =  {1, ...,  L}). 

k 

target  {k  G  1C  =  {1, ...,  iL}). 

path  of  target  k  {uk  G  Clk)- 
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Sets 

He)  Y  c 

forward  star  of  cell  c  for  searcher  of  class  1. 

Me)  c  c 

reverse  star  of  cell  c  for  searcher  of  class  1. 

V{1,  c,  d,  t) 

set  of  quadruples  {1' ,  d' ,  d" ,t')  incompatible  with  a  searcher 
of  class  1  that  moves  from  c  to  d  starting  in  time  period  t. 

Parameters 

detection  rate  in  cell  c  in  time  period  t  against  target  k  for  a 
searcher  of  class  /  when  the  searcher  previously  occupied  d . 

Cc,t{^k) 

1  if  cell  c  is  on  target  path  uik  in  time  period  f,  otherwise  0. 

number  of  searchers  of  class  1  that  occupy  cell  c  in  time 
period  0. 

Ji 

number  of  searchers  of  class  1. 

Qki^^k) 

probability  that  target  k  takes  path  ook- 

di,c,c' 

number  of  time  periods  needed  for  a  searcher  of  class  /  to 
move  directly  from  cell  c  to  cell  d  and  search  d . 

Tic 

maximum  number  of  searchers  that  occupy  cell  c  in  a  time 
period. 

Decision  Variables 

Xi^c,c',t  number  of  searchers  of  class  I  that  occupy  cell  c  in  time 

period  t  and  that  move  to  cell  d  next.  {X  denotes  the 
vector  with  components  Xi^c,c',t,  I  ^  jC,c,d  ^  C,t  &  Tq.) 

Yc^t,k  auxiliary  variable  representing  “search  effort”  in  cell  c 

against  target  k  in  time  period  t.  {Yk  denotes 
the  vector  with  components  Vc,t,fc,  c  EC,t  E  T.) 


Functions 

fkiXk)  nondetection  probability  of  target  k 


=  gfc(a;fc)exp 


^  '  Cc,t{^k)Yc^t,k 
c,t£T 


(4) 


Formulation 


min  max  fkiYk) 

(5) 

S.t.  ^  '  Xi^(J 'y  '  Xi^c,d ,t 
c'e'R.j(c)  c'^Ti{c) 

y  i,c,t  E  T 

(6) 

^  '  Xl^c,c',0  ^l,c,0 

c'eXi{c) 

y  i,c 

(7) 

^  '  'y  '  e^l,c' ,c,t,kXl^c' ,c,t—di  Y^  f  k 

l&Cc' Millie) 

y  c,t  E  T,k 

(8) 

i&Cd£ni(c) 

y  c,tET 

(9) 
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(10) 


Xl,c,c',t  +  <1  y  I,  c,c'  E  J^i{c),t, 

V  (/',  c",  c\  t')  G  V{1,  c,  c  ,  t) 

£  {0,  1,  2,  J;}  yi,C,C,t  (11) 

Yc,t,k>0  yc,tET,k.  (12) 

The  decision  variable  Fc,t,fc  could  be  eliminated  by  substitution  using  (8),  but  is  included 
for  notational  simplicity.  We  obtain  the  nondetection  probability  for  target  k  in  (4)  from  (2) 
by  application  of  the  total  probability  theorem  and  the  fact  that  detection  in  cell  c  in  time 
period  t  can  occur  only  if  the  target  occupies  that  cell  at  that  time.  The  objective  function 
(5)  aims  to  minimize  the  largest  nondetection  probability.  The  objective  function  of  SPX 
is  convex  and  the  nondetection  probabilities  fk{Yk),k  G  /C,  are  convex  and  continuously 
differentiable.  In  problem  instances  with  a  large  number  of  possible  target  paths,  calculation 
of  fkiXk)  using  (4)  is  expensive.  However,  a  Markovian  target  model  permits  a  more  efficient 
way  of  computing  fk(Yk),  which  we  present  and  exploit  in  Sections  4.2  and  4.3.  Consequently, 
SPX  is  applicable  for  both  conditional  and  Markovian  target  models. 

Constraints  (6)  and  (7)  ensure  route  continuity  and  dehne  initial  conditions  for  the 
searchers,  respectively.  Deconfliction  constraints  (9)  and  (10)  limit  the  number  of  searchers 
that  can  occupy  cell  c  to  at  most  Uc  in  any  time  period  t  E  T  and  exclude  moves  in  conflict 
with  each  other,  respectively.  We  observe  that  SPX  prescribes  the  “best”  search  plan  prior 
to  detection  of  the  first  target.  It  is  beyond  the  scope  of  the  paper  to  plan  for  events  after 
the  hrst  detection.  We  discuss  the  solution  of  SPX  in  Sections  3.3  and  4.3,  but  hrst  deal 
with  important  special  cases. 

3.2  Linearizations  of  SPX  for  Homogeneous  Searchers  and  Single  Target 

We  now  consider  special  classes  of  instances  of  SPX  involving  homogeneous  searchers  and 
a  single  target,  which  allow  us  to  construct  two  equivalent  linear  models.  For  simplicity  we 
also  exclude  deconfliction  constraints.  We  discuss  the  merits  of  extending  these  linearizations 
to  the  general  SPX  in  Section  3.3.  Specifically,  SPX  with  one  searcher  class  (i.e.,  L  =  1), 
one  target  (i.e.,  K  =  1),  a  constant  detection  rate  over  all  cells  and  time  periods,  and  no 
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deconfliction  constraints  takes  the  following  form.  In  this  case,  we  drop  the  subscripts  k  and 
I 

Model  SPl: 


Indices 

c,  c' 

cells  (c,  c'  G  C  =  {1, . . . ,  C}). 

t 

time  periods  {t  eTq  =  {0}  U  T,  T  =  {1, ...,  T}). 

u 

path  of  target  (a;  G  fl). 

Sets 

X{c)  C  C 

forward  star  of  cell  c  for  a  searcher. 

n{c)  c  c 

reverse  star  of  cell  c  for  a  searcher. 

Parameters 

a 

detection  rate  for  a  searcher. 

Cc,tA) 

1  if  cell  c  is  on  target  path  uj  in  time  period  t,  otherwise  0 

^c,0 

number  of  searchers  that  occupy  cell  c  in  time  period  0. 

J 

number  of  searchers. 

q{u) 

probability  that  target  takes  path  ui. 

number  of  time  periods  needed  for  a  searcher  to  move 
directly  from  cell  c  to  cell  d  and  search  d . 

Decision  Variables 

Xc^c',t  number  of  searchers  that  occupy  cell  c  in  time  period  t 

and  that  move  to  cell  c'  next.  {X  denotes  the  vector  with 
components  Xc,c',t,  c,c'  E  C,t  E  %.) 

Zc^t  number  of  searchers  that  occupy  cell  c  in  time  period  t. 

{Z  denotes  the  vector  with  components  Zc^t,  c  ^  C,t  &  T.) 


Functions 

f{Z)  nondetection  probability  of  target 


=  ^  q{uj)exp  -  ^  CcA^)aZc,t  ■  (13) 

Ljeo  \  c,t£T  ) 

Formulation 


min  fiZ) 

x,z  ^  ’ 

(14) 

S-t.  Xc^c'q 

c'&n(c)  ’  c'&X{c) 

y  c,teT 

(15) 

'y  Xc^c',0  ^c,0 

c'eJ^(c) 

V  c 

(16) 

^ c'  ,c,t—d^i  —  Zc,t 

y  c,teT 

(17) 

c'e7^(c) 
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Xc,c'.t  >0  V  c,c',t 
Xc.t  ^  {0»  1»  2,  \f  c,t  T 


(18) 

(19) 


The  objective  function  and  constraints  in  SPl  are  specializations  of  those  in  SPX  with 
minor  changes  caused  by  the  new  auxiliary  integer  variable  Zc,t  which  effectively  equals  Pc,i,fc 
divided  by  a.  Hence,  the  objective  function  now  involves  instead  of  Yc^t,k-  We  utilize 
the  integrality  of  Zc^t  in  the  following  linearization  as  well  as  in  construction  of  cutting  planes 
in  Section  4.1.  We  refer  to  Z  as  a  search  plan. 

We  next  derive  two  linearizations  of  SPl.  The  hrst  linearization  is  applicable  in  the 
case  of  a  conditional  target  model  with  a  moderate  number  of  possible  target  paths.  The 
second  linearization  is  limited  to  the  situation  with  a  Markovian  target  model.  We  deal  with 
the  two  linearizations  in  turn. 

The  objective  function  in  SPl  is  a  hnite  sum  of  exponential  functions  over  all  possible 
target  paths;  see  (13).  Each  exponential  function  has  as  argument  an  integer  multiple  of  a 
between  0  and  JT,  where  J  is  the  number  of  searchers.  Hence,  the  objective  function  in 
SPl  can  equivalently  be  represented  by  piecewise  linear  functions.  This  observation  leads 
to  the  hrst  linearization  of  SPl. 

Model  SPl-L; 

Indices 

As  in  SPl. 

i  number  of  looks  on  a  target  path  {i  =  0, 1, ...,  JT). 

Sets  and  Parameters 

As  in  SPl. 

Decision  Variables 

As  in  SPl. 

auxiliary  variable  representing  nondetection  probability 
given  target  path  oj. 

Formulation 

min  ^  q{uj)U^  (20) 

s.t.  e“*"(l  +  i  -  ie“")  +  -e"*“(e""  -  1)  Cc,i(^^)a^c,t  <U^  'i  u,i  (21) 

^  c,t&r 
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(15) -(19) 


SPl-L  is  a  mixed-integer  linear  program.  Constraints  (21)  ensure  that  the  optimal 
solution  results  in  a  value  of  that  is  exactly  the  conditional  nondetection  probability 
given  that  the  target  follows  path  oj.  The  other  constraints  are  identical  to  those  in  SPl. 
The  number  of  constraints  and  variables  in  SPl-L  grows  linearly  in  the  number  of  possible 
target  paths  and,  hence,  the  formulation  becomes  difficult  to  solve  for  large  number  of  such 
paths.  This  motivates  a  second  linearization  of  SPl. 

The  second  linearization  of  SPl,  denoted  SPl-LM,  assumes  a  Markovian  target  model 
where  the  target  at  time  t  G  T  moves  according  to  a  transition  probability  matrix  T*  with 
elements  7c,c',t,  c,  c'  G  C.  Specifically,  7c,c',t  is  the  probability  that  a  target  occupying  cell  c 
in  time  period  t  occupies  cell  c'  in  time  period  t  +  1.  As  we  see  below,  it  is  not  necessary  to 
enumerate  all  possible  target  paths  in  the  case  of  a  Markovian  target  model. 

We  derive  SPl-LM  from  SPl  by  introducing  an  “information  state”  Pc,t  which  equals 
the  probability  that  the  target  occupies  cell  c  in  time  period  t  and  that  the  target  has  not 
been  detected  prior  to  t.  Given  this  information  state  and  a  search  plan  with  j  searchers 
occupying  cell  c  in  time  period  f,  the  probability  of  detection  in  cell  c  in  time  period  t  and 
no  prior  detection,  is  simply 

P,,t(l-e-^“-‘),  (22) 

where  etc,*  is  the  detection  rate  of  each  searcher  in  cell  c  in  time  period  t. 

Suppose  that  a  search  plan  is  described  by  the  binary  variables  14,4 which  is  1  if  j 
searchers  occupy  cell  c  in  time  period  t  and  is  0  otherwise.  Then,  the  probability  of  detection 
over  the  full  time  horizon  becomes 


(23) 

teT  c 

The  information  state  Pc,*  depends  on  the  search  plan  as  follows.  Clearly,  =  pc,  the 
probability  that  the  target  occupies  cell  c  initially.  Moreover,  it  follows  from  the  definition 
of  Pc,t  and  the  assumption  of  a  Markovian  target  model  that 


Pc,t+i  'y  Pc'  ,c,tPc' 


(24) 
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for  all  c  G  C  and  t  =  1,  2, T  —  1.  While  the  expressions  (23)  and  (24)  are  nonlinear,  they 
can  be  linearized  as  shown  in  the  following  formulation. 

Model  SPl-LM; 


Indices 

c,  c' 
t 

j 

Sets 

As  in  SPl. 
Parameters 

^c,0 


Pc 

Qc,t 


cells  {c,  c'  E  C  =  {1, . . . ,  C}) . 

time  periods  {t  E  %  =  {0}  U  T,  T  =  {1, T}). 

number  of  searchers  in  a  cell  (j  G  JT”  =  {1, . . . ,  J}). 


detection  rate  in  cell  c  in  time  period  t  for  any  searcher, 
number  of  searchers  that  occupy  cell  c  in  time  period  0. 
number  of  time  periods  needed  for  a  searcher  to  move 
directly  from  cell  c  to  cell  c'  and  search  ch 
probability  that  a  target  that  occupies  cell  c  in  time  period  t 
occupies  cell  c'  in  time  period  t  +  1. 

probability  that  the  target  occupies  cell  c  in  time  period  1. 
probability  that  the  target  occupies  cell  c  in  time  period  t, 
i.e.,  Qc,t  X)c'  Qc' ,t—Vy c' t  2,  3, ...,  Qc,i  Pc- 


Decision  Variables 

Pc^t  probability  that  the  target  occupies  cell  c  in  time  period  t 

and  target  not  detected  prior  to  t. 

Qc,t,j  auxiliary  variable  that  equals  Pc,t(l  —  e“-^"‘=’*)  if  =  1 

and  otherwise  0. 

Vc,t,j  1  if  there  are  j  searchers  that  occupy  cell  c  in  time  period  t 

and  otherwise  0. 

Wc^t  auxiliary  variable  that  equals  if  =  1 

and  otherwise  Pc,t- 

Xc^c',t  number  of  searchers  that  occupy  cell  c  in  time  period  t  and  that 

move  to  cell  d  next. 

Formnlation 


tcET  c  j 

(25) 

d  c,t  E  T,j 

(26) 

Qc,t,j  <  (1  -  e-^“-‘)P,,i 

d  c,t  E  T,j 

(27) 

Pc,t+1  =  ^  lc',c,t^c',t 
d 

W  c,t  =  1,2,...,T-  1 

(28) 

y  c,tET 

(29) 
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c.t 


<  +  gc,i(l  -  -  K,t,,)  '^c,teT,j  (30) 

Pc,i  =Pc  V  c  (31) 

Pc,t  <  qc,t  V  c,  t  (32) 

Xc',c,t-d^,^^  =  ^jVc,t,j  yc,teT  (33) 

c'&n(c)  j 

E^c,tj<l  yc,teT  (34) 

(15),  (16) 

Pc,t>0  y  c,teT  (35) 

Qc,tj  >  0  yc,teT,j  (36) 

^  {0, 1}  y  c,t  ^  T ,j  (37) 

W,,t>0  Vc,teT  (38) 

Xc,c’,t>0  y  c,c,t  (39) 


The  objective  function  (25)  in  SPl-LM  gives  the  probability  of  nondetection;  its  cor¬ 
rectness  follows  from  (23).  However,  since  (23)  is  nonlinear,  we  linearize  it  using  the  auxiliary 
variable  Qc,t,ji  which  equals  Pc,i(l~e“-^“‘=’*)  if  Vc^t,j  =  1  and  equals  0  otherwise.  This  lineariza¬ 
tion  is  accomplished  using  constraints  (26)  and  (27).  This  is  a  “big-M”  type  of  formulation 
(see,  for  example,  [22],  pp.  642-643)  where,  in  theory,  any  constant  at  least  as  large  as  Pc^t 
would  suffice  in  front  of  (1  —  in  (26).  Recall  that  Pc^t  is  the  probability  that  the  target 

occupies  cell  c  in  time  period  t  and  target  not  detected  prior  to  t.  Moreover,  recall  that  qc,t 
is  the  probability  that  the  target  occupies  cell  c  in  time  period  t.  Hence,  gc,t  >  Pc,t  for  ah 
c,t.  Consequently,  we  set  the  “big-M”  in  (26)  to  gc,t-  We  also  use  qc,t  to  bound  the  range  of 
Pc,t  in  (32). 

The  evolntion  of  the  information  state  is  also  nonlinear;  see  (24).  In  SPl-LM,  we 
linearize  that  expression  by  means  of  the  auxiliary  variable  Wc,t  and  constraints  (28)-(30). 
Note  that  Wc^t  equals  if  Vc^t,j  =  1  and  equals  Pc^t  otherwise.  The  initial  target 

location  is  acconnted  for  in  (31).  Constraints  (33)  and  (34)  relate  Xc^c',t  to  the  binary 
variable  W.tj- 
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3.3  Linearizations  of  the  General  SPX 


SPl-L,  which  linearizes  SPl  in  the  case  of  a  conditional  target  model,  generalizes  easily 
to  a  linear  model  eqnivalent  to  SPX  if  all  detection  rates  C(i^c',c,t,k  in  SPX  are  rational 
nnmbers.  In  that  case,  all  detection  rates  SPX  can  be  expressed  as  an  integer  mnltiple  of 
a  nnmber,  say,  a.  Hence,  Yc^t,k  conld  be  expressed  as  a  times  an  anxiliary  integer  variable. 
Similar  to  the  approach  leading  to  SPl-L,  the  exponential  terms  in  fkiXk)  could  then  be 
expressed  by  piecewise-linear  functions.  Standard  techniques  for  linearizing  piecewise-linear 
functions  would  then  lead  to  a  mixed-integer  linear  program.  If  ai^c',c,t,k  differs  substantially 
across  different  elements  of  C,  C,  T,  and  1C,  a  would  need  to  be  relatively  small.  Hence, 
the  piecewise-linear  functions  may  involve  a  large  number  of  pieces  and  the  resulting  mixed- 
integer  program  may  be  large. 

Under  the  same  assumption  on  the  detection  rates  and  given  a  Markovian  target  model, 
SPl-LM  generalizes  to  a  linear  model  equivalent  to  SPX  through  a  redehnition  of  j.  While 
j  gives  the  number  of  searchers  occupying  a  cell  during  one  time  period  in  SPl-LM,  the 
new  linear  model  would  require  j  to  represent  /a,  where  denotes  the  sum  of  the 
detection  rates  of  all  searchers  occupying  a  cell  in  a  time  period.  This  sum  may  be  larger 
than  the  number  of  searchers  occupying  the  cell  as  each  searcher  would  have  a  detection 
rate  of  .^o;,  where  is  a  positive  integer.  Since  this  linearization  effectively  assigns  a  binary 
variable  to  each  possible  value  of  the  “search  effort”  in  a  cell,  the  resulting  mixed-integer 
linear  program  may  become  large. 

In  view  of  the  above  discussion,  we  see  that  linearizations  of  SPX  tend  to  be  of  rea¬ 
sonable  size  and  practical  value  when  all  detection  rates  ai^c',c,t,k  can  be  expressed  as  small 
integer  multiples  of  a.  For  example,  this  is  the  case  when  all  detection  rates  equals  1  ■  a 
for  some  a  >  0  as  in  SPl.  In  the  next  section,  we  focus  on  SPl-L  and  SPl-LM  for  this 
important  class  of  instances,  but  also  consider  the  solution  of  the  full  SPX  by  means  of  a 
cutting-plane  algorithm. 
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4  Algorithms  and  Numerical  Results 


This  section  discusses  the  solution  of  SPX,  with  particular  focus  on  SPl  as  linearizations 
appear  especially  attractive  in  that  case.  In  addition  to  linearizations,  we  also  consider  the 
cutting-plane  method  (see,  e.g.,  [15,  8,  34])  as  applied  to  SPX  and  SPl,  and  present  a 
specialized,  strengthened  cutting  plane  for  SPl.  First,  we  consider  SPl  in  the  case  of  a 
conditional  target  model  with  a  moderate  number  of  possible  target  paths.  In  practice,  this 
situation  occurs  when  the  paths  are  generated  by  Monte  Carlo  sampling  from  a  (complex) 
motion  model.  The  U.S.  Coast  Guard’s  decision  aid  GASP  generates  target  paths  in  this 
manner.  Second,  we  examine  SPl  in  the  case  of  a  Markovian  target  model.  The  Markovian 
assumption  is  adopted  in  practice,  for  example  by  the  U.S.  Navy’s  decision  aid  NODESTAR, 
as  well  as  in  the  literature  [32,  17].  We  refer  to  Chapter  7  in  [33]  for  further  details  about 
target  motion  models.  Third,  we  solve  SPX  using  a  cutting-plane  algorithm  for  the  case 
with  a  Markovian  target  model.  The  algorithm  is  also  applicable  in  the  case  of  a  conditional 
target  model,  but  we  do  not  examine  that  situation  in  detail. 

4.1  Solution  of  SPl  and  SPl-L  for  a  Conditional  Target  Model 


This  section  develops  a  specialized  cutting  plane  for  SPl  and  compares  the  performance  of 
the  resulting  cutting-plane  algorithms  with  those  corresponding  to  the  solution  of  SPl  and 
SPl-L  by  standard  solvers. 

The  standard  cutting-plane  algorithm  for  convex  (mixed-integer)  programs  sequentially 
builds  and  minimizes  successively  better  piecewise-linear  approximations  of  a  convex  func¬ 
tion;  see,  e.g.,  [15,  8,  34]  and  the  recent  review  [2].  For  completeness  and  ease  of  reference, 
we  next  state  our  implementation  of  the  standard  cutting-plane  algorithm.  Below,  we  need 
the  partial  derivatives 


df{Z) 

=  -a  g(^)w(w)exp 


(40) 


Algorithm  1  (Obtains  near-optimal  solutions  of  SPl) 


Data.  Relative  optimality  tolerances  6,6i  >  =  0, 1,  2, .... 
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Step  0.  Set  the  lower  bound,  on  the  optimal  value  of  SPl  to  0;  set  the  upper  bound, 
on  the  optimal  value  of  SPl  to  1;  and  set  i  =  1  and  =  0. 

Step  1.  Calculate  f{Z^)  and  V/(Z*).  If  /(Z*)  <  then  set  ^  =  f{Z^)- 

Step  2.  If  then  stop. 

Step  3.  Solve 

Pj  ;  min,^  (41) 

f{Zi)  +  Vf{Zl)’{Z-Z’)<i  Vj  =  l,2,...,»  (42) 

(15)  -  (19) 

to  near  optimality.  That  is,  determine  a  lower  bound  and  a  feasible  solution 
Z*+^X*+l)  of  P,  such  that  <  54*+^ 

Step  4.  If  >  4  then  set  ^ 

Step  5.  If  ^  —  ^  <  then  stop.  Else,  replace  i  by  i  +  1,  and  go  to  Step  1. 

We  note  that  Algorithm  1  is  guaranteed  to  solve  SPl  to  optimality  if  <5  =  0  and  <5*  =  0 
for  all  i.  Fixing  4  >  0  (i.e.,  accepting  near-optimal  solutions  of  Pj)  does  not  guarantee 
convergence,  but  does  improve  computational  speed.  To  balance  convergence  and  speed,  we 
allow  for  small  positive  Si  and  adopt  the  simple  safeguard  of  replacing  Si  by  (5j/2  if  a  previous 
solution  is  repeated.  This  approach  prevents  Algorithm  1  from  jamming  at  a  nonoptimal 
solution. 

The  cutting  plane  (42)  can  be  strengthened  by  taking  advantage  of  the  special  structure 
of  f(Z)  and  the  fact  that  Z^^t  is  integer.  The  strengthened  cut  uses  hnite  differences  of  the 
objective  function  f(Z)  by  considering  the  pertnrbation  from  Zc,t  to  Zc,t  +  1  while  keeping 
all  other  variables  hxed.  Theorem  1  formalizes  this  discussion,  using  Ac,t  to  denote  a  CT- 
dimensional  binary  vector  in  which  component  in  position  (c,  t)  is  1  and  the  other  components 
are  all  0. 
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Theorem  1  For  any  CT -dimensional  nonnegative  integer  vectors  Z  and  Z, 

f(Z)  +  E  ifiZ  +  A,j)  -  Jimz.,  -  4.)  <  f(Z)-  (43) 

c,t£T 

Proof.  Let  a‘^  be  a  CT-dimensional  vector  defined  by  components  Cc,i(i^)ct  and  6“^  = 
—  lng(a;).  Then,  a‘^  >0  and  6^  >  0.  Hence,  f{Z)  =  where  fuj{Z)  =  exp{—a‘^Z  — 

b‘^),  and  the  result  holds  if  f^{Z)  +  Ec,ter(/i^(^  +  ^c,t)  -  fu;{Z)){Zc,t  -  Zc,t)  <  UiZ)  for 
all  u.  Consequently,  we  need  to  show  that  /^(Z)[l  +  Ec,ter(exp(— Oct)  ~  —  Zc,t)  — 

exp(— a‘^(Z’  —  Z))]  <  0  for  an  arbitrary  target  path  ci;  G  H.  Let  f3  =  exp(— a),  and  let  M 
denote  the  set  of  the  cell-time  pairs  {c,t)  E  C  x  T  such  that  Cc,t{^)  —  ^  cell  c 

is  on  path  u  in  time  period  t).  Now,  we  only  need  to  show  that  0(/5)  =  (1  —  P)k  +  >  I, 

where  k  =  E(c,t)eAr(-^c,t  —  We  find  that  d(j){P)/dp  =  0  for  /?  =  1.  Hence,  it  follows  from 
convexity  of  0(-)  on  (0,  cx))  that  0(-)  has  a  minimum  value  of  1  for  any  k.  This  completes 
the  proof. 

We  observe  that  the  assumption  of  constant  detection  rate  a  in  SPl  is  critical  for  the 
validity  of  Theorem  1  as  the  following  counterexample  illustrates. 

Consider  a  two-cell  problem  instance  with  C  =  2,  T  =  2,  and  a  single  target  path 
w  =  (1,  2),  i.e.,  the  target  occupies  cells  one  and  two  in  time  periods  one  and  two,  respectively. 
Suppose  that  Z  =  {Zn,  Z12,  Z21,  Z22)  =  (0,  0, 1, 1)  and  Z  =  {Zn,  Z12,  Z21,  Z22)  =  (1, 1, 0,  0). 
If  the  detection  rate  is  the  same  in  cells  1  and  2,  then  it  is  easily  shown  that  equality  holds 
in  (43).  Hence,  (43)  is  satisfied  as  expected.  However,  if  the  detection  rate  is  2  in  cell  one 
and  1  in  cell  two,  then  the  right-hand  side  of  (43)  equals  exp(— 2)  ^  0.13  and  the  left-hand 
side  equals  exp(— 3)  —  exp(— 2)  -|-  exp(— 1)  ^  0.28.  Hence,  (43)  does  not  hold  in  this  case. 

We  refer  to  the  cut  of  Theorem  1  as  a  secant  cut,  and  to  Algorithm  1  with  (42)  replaced 
by 

f(zn+  Y.(fiZ‘  +  ^c.,)-f{Z’)){Zc.,-Zi,)<(  Vj  =  l,2,...,»  (44) 

c,t£T 

as  Algorithm  2. 

We  have  also  examined  the  used  of  submodular  cuts  ([19],  p.  710),  but  find  them  weak. 
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Indeed,  they  lead  to  a  cutting-plane  algorithm  that  is  10  to  100  times  slower  than  Algorithm 
2.  (New  results  in  [1]  strengthen  those  cuts  and  may  result  in  a  faster  algorithm.  We  do  not 
pursue  that  topic  here,  however.) 

We  consider  a  third  variation  of  the  cutting-plane  algorithm  where  the  continuous  re¬ 
laxation  of  Pj  is  solved  for  a  set  of  initial  iterations  after  which  Pj  is  solved.  This  version 
is  motivated  by  the  fact  that  a  large  number  of  cuts  can  be  obtained  quickly  by  solving  the 
continuous  relaxation  of  Pj.  (See  [30]  for  a  similar  idea.)  In  this  version  of  the  cutting-plane 
algorithm,  called  Algorithm  3,  (44)  is  used  when  is  integer  valued  and  (42)  is  used 
otherwise. 

We  implement  Algorithms  1-3  as  well  as  the  models  SPX,  SPl,  SPl-L,  and  SPl- 
LM  in  the  General  Algebraic  Modeling  System  (GAMS)  Distribution  22.9  [10]  on  a  laptop 
computer  with  1.0  GB  of  RAM  and  2.16  GHz  processor  running  Windows  XP.  The  mixed- 
integer  linear  programs  Pj  and  SPl-L  are  solved  (approximately)  using  GPLEX  11.2  [14] 
with  default  options. 

We  hnd  that  solution  times  of  Algorithms  1-3  are  reduced  signihcantly  when  we  only 
require  a  near-optimal  solution  of  Pj  rather  than  an  optimal  solution.  Hence,  we  normally 
run  Algorithms  1  and  2  with  (5i  =  0  and  5i  =  min{0.03,  gi/Z}  for  i  >  2,  where  gi  = 
is  computed  after  Step  1  of  iteration  i.  However,  we  use  Si  =  min{0.03,  gi/3,  5j_i/2}  if  X*  is  a 
repetition  of  a  previous  solution.  For  Algorithm  3,  we  set  (5j  =  0  if  the  continuous  relaxation 
of  Pj  is  solved  and  otherwise  follow  Algorithms  1  and  2.  The  continuous  relaxation  is  solved 
until  either  g^  <  10“^  is  achieved  or  a  user-dehned  maximum  time  is  consumed.  At  that 
point.  Algorithm  3  starts  solving  Pj.  We  use  10  minutes  as  the  maximum  time,  which  works 
well  in  our  tests. 

We  compare  Algorithms  1-3  to  the  mixed-integer  nonlinear  programming  solvers  Bonmin 
[6]  and  DIGOPT  [11]  as  implemented  in  GAMS.  DIGOPT  is  essentially  identical  to  Algorithm 
1  with  5j  =  0  for  all  i  and  an  initial  solution  of  the  continuous  relaxation  of  SPl.  DIGOPT 
uses  GPLEX  11.2  to  solve  mixed-integer  linear  programs  and  MINOS  5.51  [18]  for  nonlinear 
programs.  We  use  default  options  in  DIGOPT  with  exception  of  “stop  1,”  “maxcycles  le4,” 
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and  “epsmip  le-5,”  which  prevent  prematnre  termination. 

Bonmin  is  implemented  with  open-sonrce  solvers  Ipopt  and  Cbc  for  nonlinear  programs 
and  mixed-integer  linear  programs,  respectively,  see  [5] .  We  allow  multiple  solve  attempts  for 
Ipopt  (option  “num_retry_unsolved_random_point”  is  set  to  100),  but  otherwise  use  default 
options  for  Bonmin.  We  examine  all  available  algorithms  in  Bonmin:  BB,  OA,  QG,  Hyb, 
and  ECP.  Option  BB  is  a  branch-and-bound  algorithm  based  on  continuous  relaxation  of  the 
MINLP.  Hence,  at  each  node  in  the  branch-and-bound  tree  a  nonlinear  program  is  solved. 
(This  approach  was  examined  in  [9]  in  the  case  of  a  single  searcher.)  Option  OA  and  Ecp  are 
essentially  identical  to  DICOPT  when  applied  to  SPl,  but  with  other  solvers  for  nonlinear 
programs  and  mixed-integer  linear  programs.  Option  QG  is  an  implementation  of  the  branch- 
and-cut  algorithm  in  [21].  Option  Hyb  combines  QG  and  OA.  DIGOPT,  Algorithm  1,  and 
Bonmin,  with  OA,  QG,  Hyb,  and  EGP,  use  a  tangent  cut  of  the  form  (42),  while  Algorithm 
2  uses  the  secant  cut  (44).  The  excessive  memory  requirement  of  specialized  branch-and- 
bound  algorithms  using  bounds  based  on  expected  number  of  detections  [7,  27]  prohibits 
their  testing  on  problem  instance  of  interest  in  this  paper  [27]. 

The  calculation  times  reported  below  are  the  GPU  times  required  by  the  solvers  (“re¬ 
source  usage”  in  GAMS)  for  DIGOPT,  Bonmin,  and  GPLEX  (for  SPl-L).  For  Algorithms 
1-3,  calculation  time  is  the  total  time  required  including  cut-generation  (Step  1),  model 
generation  (Step  3),  and  mixed-integer  linear  program  solver  time  (Step  3).  Since  GAMS 
handles  cut  generation  in  Step  1  and  repeated  model  generation  (Step  3)  rather  inefficiently, 
the  fraction  of  the  total  calculation  time  used  by  GPLEX  on  Pj  is  typically  in  the  range 
0.65-0.95,  with  the  lower-end  values  dominant  for  larger  problem  instances.  Hence,  the  cal¬ 
culation  times  reported  for  Algorithms  1-3  can  be  improved,  possibly  substantially,  with  a 
more  efficient  implementation. 

We  test  algorithms  and  models  on  problem  instances  that  are  essentially  multi-searcher 
generalizations  of  those  used  in  [17,  28].  Specifically,  we  consider  a  square  AOI  with  C  =  25, 
49,  81,  121,  169,  or  225  cells.  Gells  are  numbered  from  left  to  right  and  from  top  to  bottom. 
Hence,  cell  1  is  in  the  upper-left  corner  and  cell  C  is  in  the  lower-right  corner  of  the  AOI. 
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After  each  time  period,  a  searcher  remains  in  its  current  cell  or  moves  to  a  cell  directly  above, 
below,  left,  or  right  of  the  current  cell,  if  such  a  cell  exists;  a  target  moves  similarly.  Hence, 
the  forward  and  reverse  stars  of  most  cells  consist  of  hve  cells,  except  on  the  boundary  of  the 
AOI.  For  all  allowable  cells  dec'  =  1-  The  time  horizon  T  =  7,8,  ...,  15.  These  time  horizons 
allow  the  searchers  to  have  at  least  a  moderate  chance  to  detect  the  targets.  Typically,  the 
optimal  nondetection  probabilities  are  in  the  range  0.4  to  0.8. 

We  adopt  a  conditional  target  model  and  randomly  generate  |f2|  possible  target  paths 
using  a  Markov  chain  with  a  transition  probability  matrix  dehned  as  follows.  The  probability 
of  the  target  remaining  in  a  cell  from  one  time  period  to  the  next,  denoted  p,  is  0.6,  with 
the  probability  of  moving  to  any  of  the  other  allowable  cells  being  equal.  We  use  |H|  =  1000 
unless  stated  otherwise.  The  detection  rate  of  all  searchers  is  a  =  — 3(ln0.4)/J,  where  J  is 
the  number  of  searchers  in  the  problem  instance.  This  choice  of  detection  rate  makes  one 
searcher  in  a  one-searcher  instance  as  capable  as  ten  searchers  occupying  the  same  cell  in  a 
ten-searcher  instance.  Moreover,  if  J  =  3,  then  this  detection  rate  gives  a  glimpse-detection 
probability  of  0.6.  All  searchers  occupy  cell  1  in  time  period  0  and,  hence,  can  search  cell 
1,  2,  and  \/C  -|-  1  in  time  period  1.  The  target  is  initially  located  in  the  center  of  the  AOI, 
i.e.,  Pc  =  1  if  c  =  \/C{-\/C  —  l)/2  -|-  {\/C  —  l)/2  -|-  1  and  Pc  =  0  otherwise.  These  problem 
parameters  are  similar  to  the  ones  in  [17,  28].  Testing  not  reported  here  shows  that  solution 
times  are  comparable  for  values  of  p  and  glimpse  detection  probability  in  the  range  [0.1,  0.9]. 

We  dehne  the  relative  optimality  gap  to  be  {^  —  ^)/^  for  Algorithms  1-3  and  comparable 
expressions  for  the  other  algorithms.  Table  1  shows  such  gaps  after  900  seconds  of  calculation 
time  for  instances  of  SPl  and  SPl-L  as  described  above  with  J  =  3,  C*  =  81,  and  a  varying 
number  of  time  periods.  The  time  in  seconds  to  reach  optimality  is  reported  in  brackets 
when  zero  gap  is  achieved  within  900  seconds.  When  an  algorithm  fails  to  return  a  feasible 
solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 

Bonmin  with  option  BB  fails  in  all  tests  because  no  feasible  solution  is  found;  see  Col¬ 
umn  3  of  Table  1.  These  results  indicate  that  simple  branch- and-bound  algorithms  are  not 
particular  efficient  for  SPl.  In  contrast,  cutting-plane  algorithms  (i.e.,  all  algorithms  in 
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T 

DICOPT 

BB 

OA 

QG 

Algorithms 
Hyb  Ecp 

CPLEX 

Algo.  1 

Algo.  2 

Algo.  3 

7 

0.0034 

failed 

0.0315 

0.0315 

0  [49] 

0.0315 

0[0] 

0[51] 

0  [26] 

0  [77] 

8 

0.0161 

failed 

0.0806 

0.0806 

0  [368] 

0.0806 

0[1] 

0  [253] 

0  [86] 

0  [216] 

9 

0.0247 

failed 

0.0219 

0.0221 

0.3002 

0.3002 

0[1] 

0.0147 

0  [337] 

0  [866] 

10 

0.0384 

failed 

0.4050 

0.4050 

0.4050 

0.4050 

0[12] 

0.0786 

0.0246 

0.0335 

11 

0.0592 

failed 

0.0525 

0.6838 

0.1315 

0.0600 

0  [35] 

0.1664 

0.0681 

0.0628 

12 

0.0868 

failed 

0.8555 

0.8555 

0.8555 

0.8555 

0  [174] 

0.2175 

0.1159 

0.0983 

13 

0.1052 

failed 

0.1521 

1.3821 

1.3821 

1.3821 

0.0094 

0.2968 

0.1497 

0.1425 

14 

0.1465 

failed 

1.1606 

1.1606 

0.2022 

1.1606 

0.0411 

0.4586 

0.2085 

0.2142 

15 

0.1630 

failed 

0.2933 

0.9287 

0.3463 

0.9948 

0.0650 

0.5790 

0.2711 

0.2697 

Table  1 :  Relative  optimality  gap  (^  —  0  /?  different  algorithms  after  900  seconds  of  calculation 

time  on  SPl  or  SPl-L  with  three  searchers,  81  cells,  and  a  varying  time  horizon  T.  The  time  in 
seconds  to  reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds. 
When  an  algorithm  fails  to  return  a  feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 


Table  1  with  exception  of  CPLEX  and  “BB”)  perform  better.  We  observe  that  Bonmin, 
with  cutting-plane  based  options  OA,  QG,  Hyb,  and  Ecp,  occasionally  terminates  prior  to 
900  seconds  due  to  inconsistent  performance  of  the  open-source  solvers. 


T 

DICOPT 

BB 

OA 

QG 

Algorithms 
Hyb  Ecp 

CPLEX 

Algo.  1 

Algo.  2 

Algo.  3 

7 

0.0007 

failed 

0.0347 

0.0347 

0.0347 

0.0347 

0[4] 

0  [120] 

0  [143] 

0  [137] 

8 

0.0008 

failed 

0.0654 

0.0654 

0  [782] 

0  [200] 

0[7] 

0.0003 

0  [345] 

0  [361] 

9 

0.0010 

failed 

0.0218 

0.3002 

0.0034 

0.1811 

0[12] 

0.0024 

0.0020 

0.0008 

10 

0.0037 

failed 

0.2243 

0.2243 

0.0379 

0.1053 

0  [20] 

0.0043 

0.0041 

0.0040 

11 

0.0036 

failed 

0.1028 

0.0858 

0.8304 

0.8304 

0  [531] 

0.0151 

0.0124 

0.0155 

12 

0.0049 

failed 

0.8026 

0.8026 

0.8026 

0.8026 

0.0007 

0.0356 

0.0321 

0.0343 

13 

0.0101 

failed 

1.0302 

1.0303 

0.7304 

0.7304 

0.0008 

0.0488 

0.0425 

0.0475 

14 

0.0327 

failed 

0.3311 

0.9813 

0.9813 

0.9813 

failed 

0.0920 

0.0719 

0.0755 

15 

0.0193 

failed 

0.2846 

1.1157 

0.3039 

1.0446 

failed 

0.1385 

0.1178 

0.1179 

Table  2:  Relative  optimality  gaps.  Same  problem  instances  and  parameters  as  in  Table  1  but  with 
15  searchers. 


The  solution  of  SPl-L  by  CPLEX  is  by  far  the  most  efficient  approach;  see  Column  8 
of  Table  1.  Table  1  also  illustrates  the  beneht  of  the  secant  cut  (44)  in  Algorithm  2  over  the 
tangent  cut  (42)  in  Algorithm  1.  Typically,  compared  to  Algorithm  1,  Algorithm  2  reduces 
the  solution  time,  or  the  optimality  gap  remaining  after  a  fixed  amount  of  calculation  time, 
by  a  factor  of  two.  Algorithm  3  appears  to  be  comparable  to  Algorithm  2,  but  sometimes  a 
little  faster  on  harder  problem  instances.  We  also  note  that  DICOPT  is  substantially  slower 
than  Algorithm  2  on  smaller  instances,  but  gains  the  advantage  over  Algorithm  2  as  problem 
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size  increases.  The  main  reason  for  this  effect  is  the  increase  in  overhead  associated  with  our 
implementation  of  Algorithm  2  as  problem  size  increases.  In  the  largest  example  (see  the  last 
row  of  Table  1),  DICOPT  carries  out  101  major  iterations  in  900  seconds  while  Algorithm 
2  manages  only  43  iterations.  However,  the  work  per  iteration  should  be  less  for  Algorithm 
2  (which  obtains  near-optimal  solution  of  Pj)  than  for  DICOPT  (which  obtains  optimal 
solutions  of  its  master  problems).  These  observations  illustrates  the  substantial  overhead 
associated  with  our  implementation  of  Algorithm  2  as  well  as  the  strength  of  the  secant  cut 
in  Algorithm  2  over  the  tangent  cut  used  in  DICOPT. 


|D| 

Algorithms 
DICOPT  CPLEX 

Algo.  2 

1000 

0.0161 

0[1] 

0  [86] 

2000 

0.0189 

0[10] 

0  [270] 

4000 

0.0251 

0  [48] 

0  [888] 

8000 

0.0268 

0  [177] 

0.0080 

16000 

0.0246 

0  [806] 

0.0254 

32000 

0.0269 

failed 

0.0402 

Table  3:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  (columns  2  and  4)  or  SPl-L  (column  3)  with  three  searchers,  81  cells,  time  horizon  8,  and 
varying  number  of  possible  target  paths  |D|.  The  time  in  seconds  to  reach  optimality  is  reported 
in  brackets  when  zero  gap  is  achieved  within  900  seconds.  When  an  algorithm  fails  to  return  a 
feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 


Table  2  reports  relative  optimality  gaps  for  the  same  cases  as  in  Table  1  but  with  J  =  15. 
The  advantage  of  the  secant  cut  (44)  of  Algorithm  2  as  compared  to  the  tangent  cut  (42) 
of  Algorithm  1  is  still  typically  present.  The  effect,  however,  is  much  reduced  as  Zc^t  may 
now  be  moderately  large  integers  and  hnite  differences  are  close  to  the  corresponding  partial 
derivatives.  We  find  that  the  solution  of  SPl-L  by  CPLEX  is  faster  than  the  alternatives 
as  long  as  SPl-L  is  not  too  large.  SPl-L  has  (JT  -|-  1)|D|  linearization  constraints,  which, 
after  reductions,  result  in  a  reduced  mixed-integer  linear  program  with  145,039  rows,  3,777 
columns,  and  1,142,886  nonzero  elements  for  T  =  14.  CPLEX  fails  to  hud  a  feasible  solution 
within  the  time  limit  of  900  seconds  for  this  problem  instance  and  the  one  with  T  =  15. 

We  examine  further  the  effect  of  large  instances  of  SPl-L  in  Table  3.  That  table  reports 
relative  optimality  gaps  as  the  number  of  possible  target  paths  |D|  grows.  We  see  again  that 
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c 

DICOPT 

Algorithms 
CPLEX  Algo.  2 

Algo.  3 

25 

0.2152 

0.0327 

0.2412 

0.1595 

49 

0.0927 

0  [210] 

0.1215 

0.1026 

81 

0.0396 

0[13] 

0.0246 

0.0372 

121 

0.0138 

0[3] 

0  [451] 

0.0001 

169 

0.0056 

0[1] 

0  [167] 

0  [336] 

225 

0[6] 

0[0] 

0  [106] 

0  [157] 

Table  4:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  or  SPl-L  with  three  searchers,  varying  number  of  cells,  and  time  horizon  10.  The  time  in 
seconds  to  reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds. 


C 

DICOPT 

Algorithms 
CPLEX  Algo.  2 

Algo.  3 

25 

0.5271 

0.2069 

2.8603 

0.6393 

49 

0.2505 

0.1015 

0.3954 

0.4191 

81 

0.1626 

0.0650 

0.2711 

0.3067 

121 

0.1212 

0.0191 

0.1704 

0.1820 

169 

0.0676 

0  [104] 

0.0984 

0.1292 

225 

0.0207 

0[8] 

0.0312 

0.0456 

Table  5:  Relative  optimality  gap  for  problem  instances  as  in  Table  4,  but  with  time  horizon  15. 


the  solution  of  SPl-L  by  CPLEX  is  clearly  the  fastest  approach  as  long  as  the  number  of 
possible  target  paths  is  moderate.  In  the  instances  examined,  16000  target  paths  dehne  the 
break  point.  We  also  observe  that  Algorithm  2  signihcantly  outperforms  DICOPT  unless 
the  number  of  possible  target  paths  is  large.  In  such  cases,  the  overhead  associated  with  our 
implementation  of  Algorithm  2  becomes  substantial.  For  |fl|  =  32000,  DICOPT  manages  to 
carry  out  87  major  iterations  in  900  seconds  and  reaches  an  optimality  gap  of  0.0269,  while 
Algorithm  2  completes  only  10  iterations  but  still  reaches  a  gap  of  0.0402.  This  indicates 
again  the  strength  of  the  secant  cut  used  in  Algorithm  2  as  compared  to  the  tangent  cut 
used  in  DICOPT. 

Table  4  displays  the  optimality  gaps  for  the  same  problem  instances  and  parameters  as 
in  Table  1,  but  for  a  varying  number  of  cells.  We  again  use  |D|  =  1000  possible  target  paths. 
Interestingly,  problem  instances  with  more  cells  are  easier  to  solve  than  those  with  fewer 
cells.  This  effect  results  from  the  fact  that  a  large  number  of  cells  makes  the  search  more 
difficult  as  each  possible  target  path  is  searched  only  a  relatively  small  number  of  times. 
Hence,  the  nondetection  probability  tends  to  be  large  and  the  corresponding  cuts  (42)  and 
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(44)  have  relatively  large  negative  coefficients.  This  implies  strong  cuts  in  DICOPT  and 
Algorithms  1-3,  as  well  as  relatively  few  active  constraints  in  (21).  This  effect  was  also 
eluded  to  in  [9]  for  the  case  of  a  single  searcher.  As  a  result,  problem  instances  with  225  cells 
solve  in  just  few  seconds  for  short  and  moderate  time  horizons.  Table  4  shows  that  CPLEX 
is  less  influenced  by  the  effect  than  the  algorithms  based  on  cutting  planes.  Moreover,  in 
cases  with  weak  cuts  (see  the  smaller  problem  instances  in  Table  4)  solving  the  continuous 
relaxation  of  Pj,  as  done  in  Algorithm  3,  appears  benehcial.  Table  5  show  similar  results  to 
those  of  Table  4,  but  with  a  longer  time  horizon. 


J 

DICOPT 

Algorithms 
CPLEX  Algo.  2 

Algo.  3 

1 

0.1367 

0[5] 

0  [177] 

0  [550] 

2 

0.1041 

0[13] 

0.0288 

0.0743 

3 

0.0384 

0[12] 

0.0246 

0.0335 

4 

0.0271 

0  [19] 

0.0253 

0.0277 

5 

0.0160 

0  [24] 

0.0231 

0.0203 

6 

0.0115 

0  [34] 

0.0157 

0.0149 

8 

0.0079 

0  [67] 

0.0108 

0.0105 

10 

0.0087 

0  [63] 

0.0086 

0.0084 

15 

0.0037 

0  [20] 

0.0041 

0.0040 

20 

0.0017 

0  [32] 

0.0031 

0.0040 

30 

0.0016 

0  [580] 

0.0032 

0.0042 

50 

0.0015 

0.0000 

0.0037 

0.0039 

Table  6:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  or  SPl-L  with  varying  number  of  searchers,  81  cells,  and  time  horizon  10.  The  time  in  seconds 
to  reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds. 


We  next  examine  the  effect  of  increasing  number  of  searchers.  Tables  6  and  7  report 
optimality  gaps  for  1-50  searchers  for  T  =  10  and  T  =  15,  respectively.  Again,  CPLEX 
outperforms  the  other  algorithms,  except  when  J  >  8  and  T  =  15.  Moreover,  optimality 
gaps  tend  to  decrease  as  the  number  of  searchers  increases.  This  results  from  the  fact  that 
the  continuous  relaxations  of  P*  and  SPl-L  are  stronger  for  more  searchers  because  the 
number  of  searchers  in  a  cell  may  be  a  moderately  larger  integer  in  those  cases.  (Note  that 
with  more  searchers  there  is  a  tendency  of  multiple  searchers  occupying  the  same  cell  at  the 
same  time.)  From  Tables  6  and  7  and  the  fact  that  there  are  fast  specialized  branch-and- 
bound  algorithms  available  for  single  searcher  problems  (see  [28]),  we  conclude  that  problem 
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J 

DICOPT 

Algorithms 
CPLEX  Algo.  2 

Algo.  3 

1 

0.7047 

0  [304] 

0.2283 

0.3554 

2 

0.2801 

0.0850 

0.3935 

0.4659 

3 

0.1630 

0.0650 

0.2711 

0.2697 

4 

0.0908 

0.0385 

0.1741 

0.1968 

5 

0.0712 

0.0258 

0.1396 

0.1709 

6 

0.0663 

0  0129 

0.1240 

0.1591 

8 

0.0403 

0.0092 

0.1112 

0.1408 

10 

0.0499 

failed 

0.1015 

0.1246 

15 

0.0193 

failed 

0.1178 

0.1179 

20 

0.1693 

failed 

0.1257 

0.1344 

30 

0.0093 

failed 

0.1259 

0.1251 

50 

0.7720 

failed 

0.1383 

0.1346 

Table  7:  Relative  optimality  gap  as  in  Table  6  but  with  time  horizon  15.  When  an  algorithm  fails 
to  return  a  feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 


instances  with  2-5  searchers  tend  to  be  the  most  difficult  to  solve. 

While  solving  SPl-L  using  CPLEX  may  be  the  preferred  solution  approach  in  most 
problem  instances  above,  the  approach  fails  in  cases  with  a  large  number  of  possible  target 
paths  (see  Table  3).  Hence,  we  next  consider  the  situation  with  an  extremely  large  number 
of  possible  target  paths,  but  add  the  assumption  of  a  Markovian  target  model. 

4.2  Solution  of  SPl  and  SPl-LM  for  a  Markovian  Target  Model 

This  section  empirically  compares  the  efficiency  of  solving  SPl  for  a  Markovian  target  model 
by  means  of  Algorithm  2,  with  that  of  solving  the  equivalent  linearized  model  SPl-LM  by 
means  of  CPLEX.  We  have  examined  several  cutting-plane  algorithms  for  solving  SPl,  but 
here  only  report  the  results  of  Algorithm  2,  as  they  are  typically  the  best.  We  note,  however, 
that  the  optimality  gaps  obtained  by  Algorithm  2  are  typically  half  of  those  of  Algorithm  1 
in  the  case  of  few  searchers.  Hence,  the  secant  cut  (44)  remains  superior  to  the  tangent  cut 
(42). 

A  good  cutting-plane  algorithm  requires  efficient  means  for  evaluating  f{Z),  V f{Z),  as 
well  as  the  hnite  difference  in  (44)  even  for  large  |f2|.  Since  the  expressions  (13)  and  (40)  use 
all  possible  target  paths  cu  G  12,  they  are  not  useful  in  practice  when  |f2|  is  large.  Brown  [4] 
introduces  alternative,  but  equivalent  expressions  for  f{Z)  and  V  f{Z)  that  utilize  the  fact 
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that  the  target  movement  follows  a  Markovian  target  model.  We  repeat  those  expressions 
here  with  a  slight  generalization  to  the  case  of  a  time-dependent  Markov  transition  matrix 
and  argue  how  they  lead  to  a  simple  expression  for  the  hnite  differences  used  in  (44). 

Given  a  search  plan  Z,  let  rc^t{Z)  be  the  probability  that  the  target  occupies  cell  c  in 
time  period  t  and  that  it  is  not  detected  in  time  periods  1,  2,  ...,  t  —  1,  and  let  Sc,t{Z) 
denote  the  probability  that  the  target  is  not  detected  in  time  periods  t  +  1,  t  +  2,  ...,  T 
given  that  it  occupies  cell  c  in  time  period  t.  Let  rt{Z)  =  [ri^t{Z),r2^t{Z),  ■  ■  ■  ,rc,t{Z)],  and 
let  St{Z)  =  S2,t{Z), . . . ,  sc,t{Z)].  We  dehne  rc,i(Z)  =  pc  and  Sc,t{Z)  =  1  for  any  cell 

c  G  C.  Thus,  rt{Z)  and  St{Z)  may  be  calculated  recursively  by 

rt{Z)  =  [ri,t_i(Z)exp(-aZi,j_i), . . .  ,rc,i-i(^)exp(-Q;Zc,t-i)]rt_i,  (45) 

and 

St{Z)  =  [si,t+i(Z)exp(-aZi,t+i), . . .  ,sc,i+i(^)exp(-aZc,t+i)]r(,  (46) 

where  r(  denotes  the  transpose  matrix  of  T^.  In  this  notation,  for  any  t  &  T,  we  hnd  that 

f{Z)  =  J2r,,t{Z)exp{-aZ,,t)sc,t{Z),  (47) 

cec 

and  components  of  V  f{Z)  are 

=  -arc,t{Z)  exp{-aZc,t)sc,t{Z)-  (48) 

The  calculation  of  hnite  differences  f{Z  +  ^c,t)  —  f{Z)  follows  similarly: 

f{Z  +  A,,)  -  f{Z) 

=  rc',t{^)[exp{-aZc',t  -  aA^,*)  -  exp{-aZc>,t)]sc',t{Z) 
c'ec 

=  rc,t(Z)[exp(-a(Zc,t -Fl))  -  exp(-aZc,t)]sc,i(^).  (49) 

Thus,  f{Z),  its  gradient,  and  hnite  diherence  can  be  evaluated  with  moderate  computational 
ehort.  We  utilize  (47)  and  (49)  in  Algorithm  2  in  the  following  computational  tests  on  the 
same  problems  instances  as  in  Section  4.1  and  with  the  same  computational  environment. 
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However,  we  now  consider  all  possible  target  paths  induced  by  the  Markov  chain  dehned  in 
Section  4.1. 

We  solve  SPl-LM  using  CPLEX,  with  default  options  except  for  branching  direction 
(up  hrst,  i.e.,  brdir  1)  and  the  use  of  disjunctive  cuts  (disjcuts  3),  which  was  found  to 
generally  perform  slightly  better  then  the  default  options.  Columns  2  and  3  of  Table  8  show 
relative  optimality  gaps  for  CPLEX  and  Algorithm  2  after  900  seconds  on  problem  instances 
with  three  searchers.  CPLEX  performs  signihcantly  better  when  T  <  10.  For  larger  T,  the 
big-M  formulation  of  SPl-LM  leads  to  weak  relaxations  and  long  run  times.  Columns  4 
and  5  of  Table  8  display  analogous  results  for  J  =  15.  These  instances  of  SPl-LM  become 
large,  and  CPLEX  competes  poorly  with  Algorithm  2. 


T 

J  =  3  Searchers 
CPLEX  Algo.  2 

J  =  15  Searchers 
CPLEX  Algo.  2 

7 

0[0] 

0[7] 

0.0000 

0[13] 

8 

0[1] 

0  [80] 

0.2055 

0.0002 

9 

0[7] 

0  [203] 

9.4739 

0.0021 

10 

0  [172] 

0.0287 

failed 

0.0043 

11 

0.1650 

0.0729 

failed 

0.0075 

12 

0.4964 

0.1158 

failed 

0.0106 

13 

1.1479 

0.1528 

failed 

0.0274 

14 

2.7033 

0.1737 

failed 

0.0599 

15 

21.4129 

0.2356 

failed 

0.0866 

Table  8:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  or  SPl-LM  with  three  and  15  searchers,  81  cells,  and  a  varying  time  horizon  T.  The  time  in 
seconds  to  reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds. 
When  an  algorithm  fails  to  return  a  feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 


The  effect  of  number  of  searchers,  J,  is  further  examined  in  Table  9.  As  when  |f2|  is 
moderate  (see  Tables  6  and  7),  Algorithm  2  tends  to  obtain  near-optimal  solutions  faster 
as  J  grows.  In  contrast,  SPl-LM  becomes  increasingly  harder  to  solve  as  J  increases  and 
the  resulting  model  becomes  larger.  Hence,  the  two  approaches  are  complimentary:  four 
searchers  represents  a  “cross-over  point”  above  which  Algorithm  2  has  a  clear  advantage  and 
below  which  the  solution  of  SPl-LM  prevails. 

Table  10  examines  the  effect  of  number  of  cells,  C .  As  in  Tables  4  and  5,  larger  number 
of  cells  improves  the  strength  of  the  cuts  in  Algorithm  2  and  reduces  the  optimality  gaps.  A 
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J 

Algorithms 
CPLEX  Algo.  2 

1 

0[1] 

0  [97] 

2 

0[12] 

0.0287 

3 

0  [172] 

0.0287 

4 

0.1189 

0.0320 

5 

0.3246 

0.0227 

10 

failed 

0.0074 

15 

failed 

0.0043 

Table  9:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  or  SPl-L  with  varying  number  of  searchers,  81  cells,  and  time  horizon  10.  The  time  in  seconds 
to  reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds.  When  an 
algorithm  fails  to  return  a  feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 

comparable  strengthening  of  the  continuous  relaxation  of  SPl-LM  occurs  as  the  number  of 
cells  increases.  This  allows  the  solution  of  large  problem  instances  with  a  Markovian  target 
model  in  seconds  for  small  and  moderate  T. 


C 

Algorithms 
CPLEX  Algo.  2 

25 

failed 

0.2630 

49 

0.6450 

0.1244 

81 

0  [172] 

0.0287 

121 

0[10] 

0[15] 

169 

0[1] 

0[31] 

225 

0[1] 

0  [20] 

Table  10:  Relative  optimality  gap  for  different  algorithms  after  900  seconds  of  calculation  time  on 
SPl  or  SPl-LM  with  three  searchers,  varying  cells,  and  time  horizon  10.  The  time  in  seconds  to 
reach  optimality  is  reported  in  brackets  when  zero  gap  is  achieved  within  900  seconds.  When  an 
algorithm  fails  to  return  a  feasible  solution  or  a  nonzero  lower  bound,  we  label  it  “failed.” 


4.3  Solution  of  SPX 

We  now  return  to  the  general  SPX  with  inhomogeneous  searchers,  deconfliction  constraints, 
and  multiple  targets.  As  discussed  in  Section  3.3,  it  is  straightforward  to  construct  linear 
models  of  SPX.  However,  such  models  tend  to  be  extremely  large  and  we  therefore  focus 
on  a  cutting-plane  algorithm.  We  let  “Algorithm  4”  denote  the  extension  of  Algorithm  1  to 
SPX.  Algorithm  4  is  essentially  identical  to  Algorithm  1,  except  that  the  following  problem 
is  solved  during  the  i-th  iteration  instead  of  Pp 

PX,  :  min,^  (50) 
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Time  period  0 
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□ 

Figure  1:  Area  of  interest  with  diamonds  indicating  initial  location  for  moving  targets,  asterisks 
indicating  difficult-to-search  cells,  and  circles  indicating  initial  locations  of  searchers.  The  number 
of  searchers  is  proportional  to  the  radius  of  a  circle. 

s.t.ft{Y’)  +  Vft(Y^)'{Y-Y>)<^  V  k,]  =  (51) 

(6)  -  (12). 

Moreover,  since  each  target  in  SPX  moves  independently,  fkiXk)  and  Vfk{Yk)  are  com¬ 
putable  by  extensions  of  (47)  and  (48).  We  also  extended  Algorithm  2  to  SPX  using  an 
approach  similar  to  the  one  described  in  Section  3.3.  The  resulting  algorithm  is  faster  than 
Algorithm  4  only  if  few  searchers  are  present  and  all  detection  rates  can  be  expressed  as 
small  integer  multiples  of  a  positive  constant.  Hence,  below  we  focus  on  Algorithm  4. 

We  apply  Algorithm  4  to  problem  instances  with  C  =  81  similar  to  the  instances  ex¬ 
amined  in  Sections  4. 1-4. 2;  see  Figure  1.  Cells  are  numbered  as  before.  We  consider  four 
targets  that  follow  Markovian  target  models.  At  time  period  one,  one  target  occupies  each 
of  the  cells  5,  15,  20,  and  66  (marked  with  diamonds  in  Figure  1).  After  each  time  period, 
a  target  remains  in  its  current  cell  or  move  to  a  cell  directly  above,  below,  left,  or  right  of 
the  current  cell  if  such  a  cell  exists.  The  probabilities  of  a  target  remaining  in  a  cell  from 
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one  time  period  to  the  next  is  0.4,  0.3,  0.2,  and  0.1,  respectively,  for  the  four  targets;  the 
probability  of  moving  to  any  of  the  other  allowable  cells  is  equal.  Hence,  the  target  that 
initially  occupies  cell  5  moves  slowly,  the  target  that  initially  occupies  cell  66  moves  quickly, 
and  the  other  two  targets  move  at  intermediate  speeds. 

We  consider  two  classes  of  airborne  searchers  and  set  the  travel  time  (i;^c,c'=iiiax{l, 
round ( 5c, where  round(a)  is  the  nearest  integer  to  a,  5c, c'  is  the  distance  between  c 
and  c'  measured  in  the  Euclidean  distance  between  the  center  of  the  cells,  and  z/;  is  the  speed 
of  searchers  of  class  /;  z/i  =  1,  z/2  =  2  cells  per  time  period.  Moreover,  we  let  the  forward  stars 
U  I  =  1,  2,  where  equals  the  set  consisting  of  c  and  the  four  cells 

sharing  a  side  with  c,  if  they  exist,  and  equals  the  set  of  all  cells  c'  with  di^c,c'  ^  [3,  5], 

if  they  exist.  Hence,  a  searcher  can  after  a  time  period  either  proceed  and  search  “locally” 
(i.e.,  select  a  cell  in  JE^(c))  or  transit  for  several  time  periods  to  a  distant  cell  (i.e.,  select  a 
cell  in  Tf{c)).  The  reverse  stars  7li{c),l  =  1,2,  are  dehned  similarly. 


Number  of 

Glimpse-detection  probability 

searchers 

Searcher  Class  1 

Searcher  Class  2 

Scenario 

Ji 

J2 

c  =  d 

ci^  d 

c  =  c'* 

c  A  c'* 

c  =  d 

d 

c  =  d* 

c  d  c'* 

1 

2 

1 

0.50 

0.29 

0.29 

0.16 

0.29 

0.16 

0.16 

0.09 

2 

4 

2 

0.50 

0.29 

0.29 

0.16 

0.29 

0.16 

0.16 

0.09 

3 

20 

10 

0.07 

0.03 

0.03 

0.02 

0.03 

0.02 

0.02 

0.01 

Table  11:  Scenarios  defining  problem  instances  of  SPX.  Columns  marked  with  c  =  c'  {c^  c')  give 
glimpse-detection  probability  for  a  searcher  that  occupy  (not  occupy)  the  current  cell  previously. 
An  asterisk  indicates  a  column  with  glimpse-detection  probability  for  difficult-to-search  cells;  see 
Figure  1. 


We  consider  three  scenarios  with  variable  glimpse-detection  probability  and  number  of 
searchers  as  summarized  in  Table  11.  In  scenario  1,  two  searchers  of  class  1  occupy  cell  1 
in  time  period  0  and  one  searcher  of  class  2  initially  occupies  cell  81.  The  glimpse-detection 
probability  of  a  searcher  of  the  hrst  class  is  0.50  if  the  searcher  occupied  the  current  cell  in 
the  last  time  period  (c  =  c'),  but  the  searcher’s  detection  rate  is  reduced  with  a  factor  0.5  if 
the  searcher  just  moved  into  the  cell  (c  7^  d).  In  view  of  (3),  this  implies  a  glimpse-detection 
probability  of  0.29;  see  Table  11.  This  reduction  accounts  for  the  effect,  which  we  have 
observed  in  field  experiments  with  actual  UAVs  [16,  24],  that  a  searcher  often  wastes  some 
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search  time  transiting  from  one  cell  to  another  even  if  the  cells  are  adjacent.  Using  the  model 
flexibility  of  SPX,  we  incorporate  this  effect  withont  resorting  to  a  hne  time  discretization. 

When  a  searcher  occnpies  one  of  the  cells  marked  with  an  asterisk  in  Fignre  1,  all  de¬ 
tection  rates  are  rednced  by  a  factor  of  0.5.  These  cells  represent  areas  with  poor  search 
conditions  and  conseqnently  low  detection  rates.  This  resnlts  in  a  glimpse-detection  prob¬ 
ability  of  0.29  when  c  =  d  and  0.16  when  c  7^  d .  For  the  class-2  searcher,  the  detection 
rate  is  rednced  with  a  factor  of  0.5  compared  to  class  1  in  all  sitnations,  with  resnlting 
glimpse-detection  rates  given  in  Table  11. 

Scenario  2  is  identical  to  scenario  1  except  it  has  fonr  class-1  searchers  and  two  class-2 
searchers.  Scenario  3  is  identical  to  scenario  1  except  that  there  are  20  class- 1  searchers  and 
10  class-2  searchers,  and  the  detection  rate  is  rednced  with  a  factor  of  0.1  in  all  sitnations. 
The  last  row  of  Table  11  gives  the  resnlting  glimpse-detection  probabilities.  We  note  that 
the  collective  search  effort  (i.e.,  total  detection  rate)  of  the  searchers  in  scenario  3  is  therefore 
identical  to  that  of  those  in  scenario  1.  Scenario  3,  however,  allows  more  flexibility  as  the 
search  effort  can  be  spread  more  widely. 

We  consider  both  the  sitnation  with  and  withont  deconfliction  constraints  (9)  and  (10). 
In  these  scenarios,  deconfliction  amonnts  to  ensnring  that  at  most  one  searcher  occnpies  a 
cell  each  time  period  and  that  a  searcher  is  not  allowed  to  move  from  a  cell  c  to  an  adjacent 
cell  d  G  (c)  when  another  searcher  makes  the  opposite  move  from  d  to  c.  We  assnme  that 
transit  to  a  distance  cell  (i.e.,  a  cell  d  G  T}{d))  takes  place  by  flying  at  high  altitnde,  while 
search  is  carried  ont  at  low  altitnde.  In  SPX,  we  incorporate  these  constraints  by  setting 
TT-c  =  1  for  all  c  and  'D(/,  c,  c',  f)={(/',  c",  c"',  F)  |  I'  G  £,  c"  =  c',  c"'  =  c,  F  =  t}  whenever 
d  G  jF^(c),c'  7^  c  and  otherwise  V{l,c,d,t)=(l}.  Since  searchers  transiting  between  distance 
cells  can  be  separated  easily  by  altitnde,  we  allow  the  rontes  of  snch  searchers  to  cross  each 
other  as  well  as  to  cross  over  searchers  occnpying  cells. 

Table  12  shows  lower  and  npper  bonnds  on  the  optimal  value  of  SPX  as  well  as  the 
corresponding  relative  optimality  gaps  after  15  and  60  minutes  of  calculation  time  of  Algo¬ 
rithm  4  for  these  scenarios  with  T  =  8,  10,  and  12.  Columns  4  and  5  present  the  results 
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No  deconfliction 

Deconfliction 

After 

After 

After 

After 

Scenario 

T 

Measure 

15  min 

60  min 

15  min 

60  min 

Lower  bound 

0.8514 

0.8577 

0.8519 

0.8578 

1 

8 

Upper  bound 

0.8741 

0.8663 

0.8726 

0.8663 

Relative  gap 

0.0267 

0.0100 

0.0244 

0.0099 

Lower  bound 

0.7191 

0.7267 

0.7214 

0.7310 

2 

8 

Upper  bound 

0.7521 

0.7521 

0.7736 

0.7598 

Relative  gap 

0.0458 

0.0349 

0.0724 

0.0395 

Lower  bound 

0.8367 

0.8370 

0.8875 

3 

8 

Upper  bound 

0.8374 

0.8370 

0.8875 

Relative  gap 

0.0008 

0.0000 

0  [390] 

Lower  bound 

0.7773 

0.7958 

0.7809 

0.8005 

1 

10 

Upper  bound 

0.8330 

0.8330 

0.8423 

0.8423 

Relative  gap 

0.0716 

0.0467 

0.0786 

0.0523 

Lower  bound 

0.6263 

0.6345 

0.6279 

0.6364 

2 

10 

Upper  bound 

0.7097 

0.7033 

0.7107 

0.7107 

Relative  gap 

0.1333 

0.1083 

0.1319 

0.1168 

Lower  bound 

0.7877 

0.7879 

0.8400 

0.8402 

3 

10 

Upper  bound 

0.7894 

0.7894 

0.8408 

0.8405 

Relative  gap 

0.0022 

0.0019 

0.0010 

0.0004 

Lower  bound 

0.7258 

0.7419 

0.7273 

0.7415 

1 

12 

Upper  bound 

0.8181 

0.8181 

0.8149 

0.8149 

Relative  gap 

0.1271 

0.1027 

0.1205 

0.0991 

Lower  bound 

0.5433 

0.5537 

0.5427 

0.5506 

2 

12 

Upper  bound 

0.6750 

0.6637 

0.6544 

0.6544 

Relative  gap 

0.2425 

0.1987 

0.2057 

0.1886 

Lower  bound 

0.7446 

0.7450 

0.7970 

0.7974 

3 

12 

Upper  bound 

0.7479 

0.7469 

0.7994 

0.7986 

Relative  gap 

0.0044 

0.0026 

0.0031 

0.0015 

Table  12:  Lower  and  upper  bounds  on  the  optimal  value  of  SPX  as  well  as  relative  optimality 
gaps  after  15  and  60  minutes  of  calculation  times  of  Algorithm  4  for  scenarios  1-3  with  and  without 
deconfliction  constraints.  The  time  in  seconds  to  reach  optimality  is  reported  in  brackets  when  zero 
gap  is  achieved  within  60  minutes. 
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Time  period  8 
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Figure  2:  Optimal  searcher  location  during  time  period  8  for  scenario  3,  T  =  8,  and  no  deconflic- 
tion  constraints.  The  radius  of  a  circle  is  proportional  to  the  number  of  searchers  occupying  the 
corresponding  cell  during  time  period  8  (cells  4,  7,  12,  14,  16,  20,  24,  38,  55,  66,  67,  and  75  contain 
1,  2,  1,  3,  4,  1,  3,  1,  4,  7,  2,  and  1  searchers,  respectively).  Diamonds  indicate  initial  location  for 
moving  targets  and  asterisks  indicate  difficult-to-search  cells. 

for  the  case  with  no  deconfliction  constraints,  while  columns  6  and  7  include  deconfliction 
constraints.  As  in  the  case  of  Algorithms  1-3,  Algorithm  4  solves  problem  instances  with 
more  searchers  (as  in  scenario  3)  quicker  than  those  with  fewer  searchers  (as  in  scenario  1). 
We  also  find  longer  time  horizons  to  be  more  difficult,  again  primarily  due  to  the  weaker 
cuts  in  the  case  of  smaller  nondetection  probabilities. 

Deconfliction  constraints  restrict  SPX  and  result  in  an  increase  in  the  optimal  value. 
In  scenarios  1  and  2,  the  change  is  small  due  to  the  relatively  low  number  of  searchers. 
Deconfliction  constraints  increase  the  optimal  value  with  about  0.05  in  scenario  3  where  30 
searchers  are  present.  Hence,  deconfliction  constraints  effectively  force  the  searchers  to  give 
up  0.05  in  probability  of  detection. 

Figure  2  illustrates  the  optimal  location  of  searchers  in  time  period  8  for  scenario  3, 
T  =  8,  and  no  deconfliction  constraints.  The  radius  of  a  circle  is  proportional  to  the  number 
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Time  period  8 
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Figure  3:  Optimal  searcher  location  during  time  period  8  for  scenario  3,  T  =  8,  and  deconfliction 
constraints.  A  dot  indicates  one  searcher. 

of  searchers  located  in  the  corresponding  cell  during  time  period  8  (cells  4,  7,  12,  14,  16,  20, 
24,  38,  55,  66,  67,  and  75  contain  1,  2,  1,  3,  4,  1,  3,  1,  4,  7,  2,  and  1  searchers,  respectively). 
We  observe  that  multiple  searchers  focus  on  a  relatively  small  number  of  cells  with  high 
probability  of  containing  targets.  This  tendency  is  also  present  in  earlier  time  periods.  In 
fact,  the  ten  searchers  initially  occupying  the  bottom-right  corner  of  the  AOI  (see  Figure 
1)  ignore  cells  near  their  initial  location,  make  a  long-distance  move  to  the  nearest  target 
(cell  66),  and,  due  to  their  high  speed,  arrive  there  for  search  in  time  period  3.  Seven 
of  the  ten  searchers  remain  in  cell  66  for  the  remaining  time  periods  to  avoid  the  lower 
detection  rate  associated  with  moving  to  another  cell.  Three  searchers  fan  out  to  adjacent 
cells,  but  noticeably  avoid  the  difficult-to-search  cells  marked  with  asterisks  in  Figure  2.  The 
20  searchers  starting  in  the  upper-left  corner  (see  Figure  1)  spread  more  widely,  even  tackle 
difficult-to-search  cells,  but  remain  somewhat  clustered. 

Figure  3  shows  the  optimal  location  of  searchers  in  time  period  8  as  in  Figure  2,  but  in 
the  case  with  deconfliction  constraints.  Now,  each  cell  allows  at  most  one  searcher  and  the 
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searchers  are  forced  to  spread  out.  As  seen  from  row  10  in  Table  12,  spreading  the  effort  is 
less  effective  and  the  optimal  value  increases  from  0.8370  to  0.8875. 


T 

No  deconfliction 
UB  Rel.  gap 

Deconfliction 

UB  Rel.  gap 

10 

0.7894 

0.0019 

0.8405 

0.0004 

14 

0.7093 

0.0038 

0.7600 

0.0016 

18 

0.6433 

0.0063 

0.6906 

0.0023 

22 

0.5869 

0.0117 

0.6292 

0.0021 

26 

0.5393 

0.0225 

0.5755 

0.0041 

30 

0.5126 

0.0865 

0.5266 

0.0053 

34 

0.4755 

0.1171 

0.4833 

0.0093 

38 

0.4473 

0.1811 

0.4452 

0.0222 

42 

0.4102 

0.2337 

0.4114 

0.0598 

Table  13:  Upper  bounds  (UB)  on  the  optimal  value  of  SPX  as  well  as  relative  optimality  gaps 
after  60  minutes  of  calculation  times  of  Algorithm  4  for  scenario  3  with  varying  T,  and  with  and 
without  deconfliction  constraints. 


Table  13  further  examines  the  upper  bounds  (UB)  on  the  optimal  value  of  SPX  and 
corresponding  relative  optimality  gaps  in  scenario  3  as  T  increases.  As  in  Table  12,  we  hnd  a 
worsening  in  solution  quality  after  60  minutes  of  calculation  time  as  T  increases.  However,  the 
increase  is  moderate  and  essentially  insignihcant  for  the  case  with  deconfliction  constraints. 
We  are  able  to  obtain  near-optimal  solutions  even  for  long  time  horizons.  Interestingly, 
deconfliction  constraints  reduce  optimality  gaps  in  these  instances  even  though  they  increase 
the  model  size. 


5  Conclusions 


This  paper  has  presented  models  and  algorithms  for  a  discrete-time  route-optimization  prob¬ 
lem,  denoted  SP,  where  multiple  searchers  and  one  or  more  probabilistically  moving  targets 
operate  on  a  hnite  set  of  cells.  We  have  formulated  a  novel  convex  mixed-integer  nonlinear 
program  (MINLP)  for  SP  that  considers  searcher  deconfliction,  time  of  search,  and  target- 
and  location-dependent  search  effectiveness.  MINLP  allows  modeling  of  real-world  situations 
not  previously  considered  including  those  with  many  searchers. 

We  propose  two  solution  approaches  for  the  MINLP  that  result  in  the  hrst  practical, 
exact  algorithms  for  SP.  One  approach  is  based  on  the  cutting-plane  method  and  leads  to 
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near-optimal  solutions  with  5%  relative  optimality  gap  in  less  than  15  minutes  of  calculation 
time  for  problem  instances  with  three  searchers,  81  cells,  and  10  time  periods.  We  enhance 
the  approach  by  developing  a  specialized  secant  cut  that  reduces  the  solution  time  with  a 
factor  of  two  in  direct  comparisons  with  a  standard  tangent  cut  on  difficult  instances  of  the 
MINLP.  In  all  problem  instances  examined,  existing  branch-and-bound  algorithms  either  fail 
to  find  a  feasible  solution  within  15  minutes  or  require  memory  in  excess  of  one  gigabytes. 

The  other  approach  is  based  on  two  novel  linearizations  of  the  MINLP  available  in 
important  special  cases  involving  a  single  target.  In  the  case  of  a  moderate  number  of 
possible  target  paths,  a  linearization  is  easily  solved  by  standard  solvers  when  the  number 
of  searchers  is  less  than  10.  This  leads  to  reduction  in  solution  times  with  one  or  two  orders 
of  magnitude  when  compared  to  the  cutting-plane  approach  and,  for  example,  allows  the 
solution  of  a  problem  instance  with  three  searchers,  81  cells,  and  10  time  periods  in  12 
seconds.  In  the  case  of  a  target  moving  according  to  a  Markov  chain,  the  solution  times  for 
solving  a  linearization  of  MINLP  is  longer  but  an  instance  with  three  searchers,  81  cells,  and 
10  time  periods  is  still  solved  to  optimality  in  less  than  three  minutes. 

The  cutting-plane  approach,  and  to  a  less  extent  the  linearization  approach,  scale  well 
as  the  number  of  searchers  grows.  Empirically,  we  observe  that  the  cutting-plane  approach 
exhibits  an  approximately  constant  solution  time  to  reach  near-optimality  as  the  number 
of  searchers  grows.  This  enables  the  computation  of  near-optimal  solutions  of  problem 
instances  with  up  to  50  searchers  in  less  than  15  minutes.  In  a  realistic  scenario  involving 
30  searchers  divided  in  two  classes  with  different  speed  and  sensor  quality,  four  targets  of 
variable  characteristics,  deconfliction  constraints,  81  inhomogeneous  cells,  and  a  long  time 
horizon  of  38  periods,  we  obtain  a  3%  near-optimal  solution  in  30  minutes.  Hence,  we 
overcome,  in  part,  the  difficulty  of  solving  problem  instances  with  many  searchers. 

Problem  instances  remain  challenging  when  they  have  two  or  more  classes  of  searchers, 
each  with  only  a  few  searchers,  and  when  long  time  horizons  lead  to  low  nondetection  prob¬ 
abilities.  A  heuristic  allocation  of  searchers  to  nonoverlapping  subareas  may  be  a  practical 
approach  in  the  first  situation  and  a  rolling-time  horizon  and/or  fix-and-relax  optimization 
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may  heuristically  address  the  second  situation.  However,  fast  solution  of  these  classes  of 
instances  of  SP  remains  a  goal  for  future  research. 
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